Critical behaviour of a surface reaction model 
with infinitely many absorbing states 



Iwan Jensen 

^ ■ Department of Mathematics, The University of Melbourne, 

; Parkville, Victoria 3052, Australia. 



Q 



e-mail: iwan@mundoe.maths.mu.oz.au 
February 1, 2008 



> 
in 
\o 
o 

CN ■ Abstract 

m 

o\ ■ 

In a recent letter [J. Phys. A 26, L801 (1993)], Yaldram et al studied 
\ the critical behaviour of a simple lattice gas model of the CO-NO catalytic 

reaction. The model exhibits a second order nonequilibrium phase transi- 
I tion from an active state into one out of infinitely many absorbing states. 

Q ' Estimates for the critical exponent [3 suggested that the model belongs to 

O ■ a new universality class. The results reported in this article contradict 

this notion, as estimates for various critical exponents show that the model 
^ . belongs to the universality class of directed percolation. 

PACS Numbers: 05.70.Ln, 05.50. +q, 64.90. +b 
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Nonequilibrium phase transitions occur in many models studied in physics, chem- 
istry, biology or even sociology. A special group of models, that have attracted a 
great deal of interest in recent years, exhibit a continuous transition into an ab- 
sorbing state. The best known examples are probably directed percolation (DP) 
Reggeon field theory [|], ||, the contact process 0-@], and Schlogl's first 
and second models Extensive studies of these and many other models 



T^-|2^ with a unique absorbing state have revealed that they belong to the 



same universality class. This provides firm support for the conjecture that con- 
tinuous transitions into a unique absorbing state generically belong to the DP 



class fO, O 



For models with multiple absorbing states the situation is not so simple. Some 
studies of two-dimensional surface reaction models yield critical exponents differ- 
ent from those of directed percolation in (2-1-1 )-dimensions p3| , p^ . However, the 
pair contact process (PCP) and dimer reaction (DR) model (in one dimension) 
clearly belongs to the DP universality class at least as far as the static 

critical behaviour is concerned. In all of these models the number of absorbing 
configurations grows exponentially with system size. However, all of the absorb- 
ing configurations are characterized by the vanishing of a unique quantity, e.g., 
the number of particle pairs in the PCP or in other cases the p3[ ^ number of 
nearest neighbor vacancy pairs. 

Recently, Yaldram et al. ||2^, studied the critical behaviour of a simple lattice 
model of the CO-NO catalytic reaction in which CO NO ^ CO2 + ^Na. 
Schematically the reaction steps are given as: 



C0^ + * 
NO^ + 2* 
CO" + O" 



CO'', 

col + 2*, 
+ 2*, 



(1) 
(2) 
(3) 
(4) 



where the superscript g (a) refers to a molecule in the the gas fase (adsorbed on 
the surface) and * marks an empty site. The catalytic surface is modelled by 
a two-dimensional triangular lattice. The rules of the computer algorithm are 
quite simple, with probability p a CO molecule is adsorbed on an empty site and 
with probability 1 — p NO adsorption is attempted. Since NO dissociates upon 
adsorption it requires a nearest neighbour pair of empty sites. In the simulations 
this is done by first chosing an empty site at random and then chosing one of the 
six nearest neighbours randomly, if the neighbour is empty O is placed on the 
original site and N on its neighbour. After each adsorption the nearest neighbours 
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are checked (in random order) and CO+0 reacts to form CO2 which leaves the 
surface at once, hkewise N+N forms N2 which desorbs immediately. It is thus 
obvious that any state without empty sites is absorbing. All processes depend 
on the presence of empty sites so an efficient algorithm uses a list of these. After 
each attempted adsorption the time variable is incremented by 1/Nf,, where 
is the number of empty sites prior to the attempt, thus making each time step 
equal to (on the average) one attempted update per lattice site. The algorithm 
outlined above differs from that used by Yaldram et al. in one aspect, when NO 
adsorbs they choose a pair of empty sites at random, whereas I choose one empty 
sites and a nearest neighbor and only adsorb NO if the nearest neighbor is empty. 
This makes NO adsorption less likely in my algoritm. However, one would expect 
this merely to lead to a change in the location of the phase transitions not to a 
change in the critical behaviour. Computer simulations by Yaldram et al. [27| 



shows that when p < pi the system always enters an absorbing state in which the 
lattice is covered by a mixture of O and N (but off course without any nearest 
neigbour pairs of N). Note that the symmetry of the lattice prevents a CO from 
being surrounded by N, as some of these N would have to be nearest neighbours 
and thus react. The number of absorbing configurations grows exponentially 
with system size. Note also that an absorbing configuration, though not unique, 
is characterized by the vanishing of the number of empty sites. At pi the model 
exhibits a continuous phase transition into an active state in which the catalytic 
process can prodeed indefinetely. Finally when p exceeds a second critical value 
P2 the model exhibits a discontinuous phase transition into a CO and N covered 
state. The phase diagram of the CO-NO reaction model is thus very similar to 



that found in various similar catalytic model |Tj, |2^, Near the critical point 
Pi one would expect the concentrations px of various lattice sites X {X = O, N, 
CO, or an empty site) to follow simple power laws. 



px-pf cx(p-pi)^-, (5) 

where is the saturation concentration. Note that the saturation concentration 
for empty sites and CO is zero, whereas it is non-zero for O and N. Yaldram et al. 
27[] found that pi = 0.185(2), where the figure in parenthesis is the uncertainty in 



the last digit, and jSx = 0.20 — 0.22. The estimates for j3x are much smaller that 
the value f3 = 0.592(10), obtained using the scaling relation (3 = 6i>\\ |TT] with 



6 = 0.460(6) and z/y = 1.286(5) [^, for directed percolation in (2+l)-dimensions. 
This could indicate that the CO-NO model belongs to a new universality class. 
However, the uncertainty in the estimate for pi is quite large, especially con- 
sidering that the /3 estimates are obtained using values of p — pi between 0.01 
and 0.001, which overlaps the error estimate for pi. Moreover, the lattice sizes 
(40 X 40) used in the simulations are very small. Actually for such small lattice 
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sizes one would expect finite-size effects to be quite prominent. All in all I think 
there is ample reason to doubt the validity of the exponent estimates obtained 
by Yaldram et al. 



In this article I report the results of extensive simulations of the CO-NO model 
using time-dependent simulations and finite-size scaling. The general idea of 
time-dependent simulations is to start from a configuration which is very close to 
the absorbing state, and then follow the "average" time evolution of this config- 
uration by simulating a large ensemble of independent realisations. This method 
is straight forward and very successful for models with a unique absorbing state 



TT| , p!5| , p!8| , p!9| , |28[| . For models with multiple absorbing state the situation is 
more intricate, as a recent study revealed that the dynamic critical exponents 
predicted via time- dependent simulations depend upon the choice of initial con- 
figuration. However, two important facts emerged from this study, first of all the 
predictions for the location of the critical point was allways correct, and secondly 
if one uses an initial configuration reminiscent of a typical absorbing configuration 
the predictions for the dynamical critical exponents coincide with those expected 
from the static critical behaviour. A recent more thorough study by Mendes et al. 

29(1 have confirmed this picture and led to a generalized scaling ansatz for models 



with multiple absorbing states. In this study I generate the initial configuration 
by simulating the CO-NO model on a 128x128 lattice (with periodic boundary 
conditions) at the value of p under investigation until it enters an absorbing state. 
An off-set (x, y) is then chosen randomly on this lattice. Hereafter the configura- 
tion is mapped cychcally onto a larger (512x512) lattice such that (x, y) is at the 
origin of the larger lattice. The particle at position (i, j) on the large lattice is the 
same as the particle at position {i+x mod 128, mod 128) on the small lattice. 
Hereafter a pair of empty sites is placed at the origin. The size of the large lattice 
ensures that the cluster of empty sites grown from the seed at the origin never 
reaches the boundaries of the lattice. We thus start in a configuration close to an 
absorbing state (just two sites are open) and it should be close to a typical ab- 
sorbing state of the infinite system. For each such configuration I simulated 5000 
independent samples and typically 50-100 independent configurations for a total 
of 250-500,000 samples. Each run had a maximal duration of 2000 time steps, 
but most samples enters an absorbing state before this limit is reached. As usual 
in this type of simulation I measured the survival probability -P(t), the average 
number of empty sites n{t), and the average mean square distance of spreading 
-R^(t) from the origin. Notice that n{t) is averaged over all runs whereas R^{t) 
is averaged only over the surviving runs. In accordance with the scaling ansatz 
for models with a unique absorbing state [jll], ^ it follows that these quantities 
have the following scaling form. 
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Pit) oc r^$(At^/^ii), (6) 
nit) oc t''^(Ati/'^ii), (7) 
R\t) cx t^QiAt^/^w), (8) 



where A = — pi| is the distance from the critical point, and i/y is the time-hke 
correlation length exponent. If the scaling functions $, \E', and G are non-singular 
at the origin it follows that Pit), nit), and -R^(t) behave as power-laws at pi with 
critical exponents —6, t], and z, respectively, for t — > oo. Generally one has 
to expect corrections to a pure power law behaviour so that, e.g., P(t) is more 
accurately given by ||2^ 



Pit) oc t-\l + at'^ + ht-^' + ■■■) (9) 

and similarly for nit) and P?it). More precise estimates for the critical exponents 
can be obtained if one looks at local slopes 



log[P(t)/P(t/m)] 

-oit) = — , (10) 

log(m) 

and similarly for r]it) and z(t). In a plot of the local slopes vs 1/t the critical 
exponents are given by the intercept of the curve for pi with the y-axis. The off- 
critical curves often have very notable curvature, i.e., one will see the curves for 
p < Pi veering downward while the curves for p > pi veer upward. This enables 
one to obtain accurate estimates for pi and the critical exponents. In Fig. 1 I 
have plotted the local slopes for various values of p. From the plot of rjit) it is 
clear that the two lower curves, corresponding to p = 0.1781, and 0.1782, veers 
downward showing that pi > 0.1782. Likevise the upper curve, p = 0.1785, has a 
pronounced upwards curvature. Though it is less evident it also seems that the 
curve for p = 0.1784 veers upwards. All in all I conclude that pi = 0.1783(1). 
This estimate differs quite a bit from that of Yaldram et al. ipi = 0.185(2)), 
which is probably due to the slightly different algorithms. Note that since NO 
adsorption is less efficient in my algorithm one would expect my estimate for pi to 
be smaller, as is also observed in the simulations. From the intercept of the critical 
curves with the y-axis I estimate 6 = 0.45(1), rj = 0.220(5) and z = 1.12(1). 
These values agrees very well with those obtained from computer simulations of 



directed percolation in (2+l)-dimensions [28|, 6 = 0.460(6), rj = 0.214(8) and 
z = 1.134(4). 
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From these results it seems reasonable to conclude that the CO-NO model belongs 
to the DP universality class. However, due to the somewhat arbitrary choice of 
the initial configuration employed in the time-dependent simulations it would be 
nice to validate this conclusion through other means. To this end I have also 
performed extensive steady-state simulations using a finite-size scaling analysis. 
Finite-size scaling, though originally developed for equilibrium systems, is also 
applicable to nonequilibrium second-order phase transitions as demonstrated by 
Aukrust et. al. [0. Their method was later applied to models with infinitely 
many absorbing states [^, As in equilibrium second-order phase transitions 
one assumes that the (infinite-size) nonequilibrium system features a length scale 
which diverges at criticality as, ^(p) oc A"'"-^, where is the correlation length 
exponent in the space direction. The basic finite-size scaling ansatz is that the 
various quantities depend on system-size only through the scaled length L/^, or 
equivalently through the variable AL^/'^^, where L is the linear extension of the 
system. Thus we assume that the density of empty sites (which will be used as 
the order parameter of the model) depends on system size and distance from the 
critical point as: 

Ps{p,L)o,L-^'^^T{^L^/^^), (11) 

such that at pi 

Ps{puL)<xL~P/''\ (12) 

In ps, and other quantities, the subscript s indicates an average taken over the 
surviving samples. Fig. 2 shows a plot of the average concentration of particles 
log2[ps(pi, L)] as a function of loggL at the critical point, pi = 0.1783. All 
simulations were performed on lattices of size L x L using periodic boundary 
conditions. The maximal number of timesteps in each trial, Im, and number 
independent samples, Ns, varied from = 300, Ns = 50,000 for L = 8 to 
tM = 125, 000, Ns = 500 for L = 256. The slope of the line drawn in the figure 
is /?/z/_L = 0.81, which comes from the DP estimate f3/i'± = 0.81(2), using the 
earlier cited estimate for f3 and z/^ = 0.729(8) p8|. The data falls very nicely on 
the line drawn using the DP estimate thus confirming that the model belongs to 
the DP universality class. 

Near the critical point the order parameter fiuctuations grow like a power law, 
Xs = L'^iip'^) — (p)^) oc A''', from which we expect the following finite-size scaling 
form, 

X.(p,L)ocL^/^-^(ALi/^-), (13) 
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such that at pi 

Xs{pi.L)^L^I''K (14) 

Fig. 3 shows a plot of \og2[Xs{pi^ L)] vs logj The slope of the straigth line 
is 0.39 as obtained from the DP value 7/z/_l = 0.39(2), where I used that 7 = 
— = 0.285(11) with 7^^ = 1.571(6) [^. The excellent agreement between 
the data and the DP-expectation confirms the DP critical behaviour of this model. 

One expects a characteristic time for the system, say the relaxation time, to scale 
like 



t{p, L) oc L-'^ii/"^T(AL1/'^^), (15) 

such that at pi 

r(pi,L) oc L^^'ii/^^. (16) 

In Fig. 4 I have plotted log2[r/j(pi, L)], where Th is the time it takes for half 
the samples to enter an absorbing state, as a function of log2 L. The slope of 
the line drawn in the figure is vw/vi. = 1.764, as obtained from the DP estimate 



28|| i'\\/iy± = 1.764(7). The DP estimate is derived from the scaling relation 
= 2/z using the earlier cited estimate for z. As can be seen the data for 
the CO-NO model is again fully compatible with DP critical behaviour. 

One may also study the dynamical behaviour by looking at the time dependence 
of Ps{pi, L, t). For t^l and L ^ 1 one can assume a scaling form 



p,(pi,L,t) oc L'^'^^Hit/U^^/''^). (17) 

At pi the system shows a power law behaviour for t < before finite-size 

effects become important. Thus for L ^ 1 and t < L'^n/'^^, p(pi, L, t) oc t^^. From 
Eq. (0) we see that this is the case for large L only if ^ = /^/z/y . It can be shown 
ITTH that this ratio also equals the critical exponent 6. Fig. 5 shows the short-time 
evolution of the concentration of empty sites at pi with L = 256, tM = 10,000, 
and Ns = 1000. The asymptotic behaviour is consistent with 6 = 0.45, as seen 
from the slope of the line. This estimate agrees well with the value for directed 
percolation 6 = 6 = 0.460(6), or the estimate 6 = 0.45(1) obtained from the 
time-dependent simulations presented above. 

In conclusion, we have provided very convincing evidence that the critical expo- 
nents of the two dimensional CO-NO model are the same as those of directed 
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percolation in (2+l)-dimensions. This is the first time that a two dimensional 
multi-component model with infinitely many absorbing states has been firmly 
placed in the DP universality class. This results lends further support to the 
extensions of the DP conjecture to models with multiple components |Q and/or 
infinitely many absorbing states at least in cases where the absorbing 



states can be characterized by the vanishing of a unique quantity. 
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Figure Captions 



Figure 1 Local slopes —S(t) (upper panel), T]{t) (middle panel), and z{t) (lower 
panel), as defined in Eq. |TU] with m = 5. Each panel contains five curves with, 
from bottom to top, p = 0.1781, 0.1782, 0.1783, 0.1784 and 0.1785. 

Figure 2 The concentration of empty sites \og2[ps{pi, L)] vs log2 L. The slope of 
the straight line is = 0.81. 

Figure 3 The fiuctuations in the concentration of empty sites \og2[Xs{Pi, L)] vs 
logg L. The slope of the straight line is ■j/i'± = 0.39. 

Figure 4 The time before half the samples enter an absorbing state log2[Th{pi, L)] 
vs \0g2L. The slope of the straight line is i^\\/i^± = 1.764. 

Figure 5 Log-log plot of ps{p,L,t), for p = pi = 0.1783 and L = 256, as a 
function of t. The slope of the straight line is ^ = 0.45. 
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